clear; close all;

saveArg = 1;
listFile = dir('Exp_*');
Lb = ([77/10, 79/10, 85/100, 60/100, 69/10]/9.8).^(2/3);


%% Reading data

accrueL = nan(40,length(listFile));
accrueB = accrueL;
accrueTime = accrueL;
accrueErupt = accrueL;
accrueThetaMean = accrueL;
accrueThetaDev = accrueL;
accrueThetaMedian = accrueL;
accrueThetaMad = accrueL;
accrueCorrectedMean = accrueL;
accrueCorrectedDev = accrueL;
accrueCorrectedMedian = accrueL;
accrueCorrectedMad = accrueL;
accrueReducedMean = accrueL;
accrueReducedDev = accrueL;
accrueReducedMedian = accrueL;
accrueReducedMad = accrueL;
accrueReducedCorMean = accrueL;
accrueReducedCorDev = accrueL;
accrueReducedCorMedian = accrueL;
accrueReducedCorMad = accrueL;
accrueVelocity = accrueL;

for i = 1:length(listFile)
    load(sprintf('%s\\Output.mat',listFile(i).name));
    
    accrueL(1:length(L),i) = L;
    accrueB(1:length(L),i) = B;
    accrueTime(1:length(L),i) = time;
    accrueErupt(1:length(L),i) = intrudeErupt;
    accrueVelocity(1:length(L),i) = velocityMean;
    
    accrueThetaMean(1:length(L),i) = thetaMean;
    accrueThetaDev(1:length(L),i) = thetaDev;
    accrueThetaMedian(1:length(L),i) = thetaMedian;
    accrueThetaMad(1:length(L),i) = thetaMad;
    
    accrueCorrectedMean(1:length(L),i) = correctedMean;
    accrueCorrectedDev(1:length(L),i) = correctedDev;
    accrueCorrectedMedian(1:length(L),i) = correctedMedian;
    accrueCorrectedMad(1:length(L),i) = correctedMad;
    
    accrueReducedMean(1:length(L),i) = reducedMean;
    accrueReducedDev(1:length(L),i) = reducedDev;
    accrueReducedMedian(1:length(L),i) = reducedMedian;
    accrueReducedMad(1:length(L),i) = reducedMad;
    
    accrueReducedCorMean(1:length(L),i) = reducedCorMean;
    accrueReducedCorDev(1:length(L),i) = reducedCorDev;
    accrueReducedCorMedian(1:length(L),i) = reducedCorMedian;
    accrueReducedCorMad(1:length(L),i) = reducedCorMad;
end

clear reduced* corrected* theta* vector* B L intrudeErupt time ...
    x y vx vy velocityMean indDir


%% Displaying

figure(1); clf;
set(gcf,'Units','normalized','OuterPosition',[0 0.05 1 0.95]);

subplot(2,2,1);
data = importdata('Exp_Prev_C1\Data_03-1800-1900\B00001.dat');
data = data.data;
vx = data(:,3);         vy = data(:,4);
noInd = find(vx==0&vy==0);
vx(noInd) = [];         vy(noInd) = [];
vector = sqrt(vx.^2+vy.^2);
theta = acos(vx./vector);
hist(theta,100); hold on;
ylim([0,70]);
xlim([0,pi]); set(gca,'XTick',0:pi/2:pi,'XTickLabel',[]);
plot(accrueReducedMedian(3,3)*[1,1],get(gca,'YLim'),'k-','LineWidth',2);
plot((accrueReducedMedian(3,3)-accrueReducedMad(3,3))*[1,1],...
    get(gca,'YLim'),'k--','LineWidth',2);
plot((accrueReducedMedian(3,3)+accrueReducedMad(3,3))*[1,1],...
    get(gca,'YLim'),'k--','LineWidth',2);
tx = text(10/180*pi,62,'a'); set(tx,'FontSize',30);
plot([1.6987,2.3902],50*[1,1],'k-','LineWidth',2);
plot([2.3,2.3902],[49,50],'k-','LineWidth',2);
plot([2.3,2.3902],[51,50],'k-','LineWidth',2);
tx = text(120/180*pi,60,'\sigma'); set(tx,'FontSize',30);
set(gca,'Fontsize',18);
set(gca,'Position',[0.08 0.59 0.38 0.36]);

subplot(2,2,3);
data = importdata('Exp_Prev_C1\Data_22-9400-9500\B00001.dat');
data = data.data;
vx = data(:,3);         vy = data(:,4);
noInd = find(vx==0&vy==0);
vx(noInd) = [];         vy(noInd) = [];
vector = sqrt(vx.^2+vy.^2);
theta = acos(vx./vector);
hist(theta,100); hold on;
ylim([0,1400]);
xlim([0,pi]); set(gca,'XTick',0:pi/2:pi,'XTickLabel',0:90:180);
plot(accrueReducedMedian(22,3)*[1,1],get(gca,'YLim'),'k-','LineWidth',2);
plot((accrueReducedMedian(22,3)-accrueReducedMad(22,3))*[1,1],...
    get(gca,'YLim'),'k--','LineWidth',2);
plot((accrueReducedMedian(22,3)+accrueReducedMad(22,3))*[1,1],...
    get(gca,'YLim'),'k--','LineWidth',2);
xlabel('\theta (deg)');
tx = text(10/180*pi,1200,'b'); set(tx,'FontSize',30);
set(gca,'Fontsize',18);
set(gca,'Position',[0.08 0.15 0.38 0.36]);



subplot(2,2,[2,4]);

xr = accrueL/1000./repmat(Lb,size(accrueL,1),1);
yr = accrueReducedMad;

for i = 1:length(listFile)
    plot(xr(:,i),yr(:,i)*180/pi,'-',...
        'LineWidth',2,'Color',[1-i/5,1-i/5,1]); hold on;
end
xlim([0.1,4]); ylim([9,60]);
plot(get(gca,'XLim'),35*[1,1],'k--');
tx = text(1.4,37,'circulating'); set(tx,'Fontsize',16);
plot(2.1*[1,1],[39,42],'k','LineWidth',1);
plot(2.1*[0.99,1],[41,42],'k','LineWidth',1);
plot(2.1*[1.01,1],[41,42],'k','LineWidth',1);
tx = text(1.05,33,'buoyant/erupting'); set(tx,'Fontsize',16);
plot(2.1*[1,1],[31,28],'k','LineWidth',1);
plot(2.1*[0.99,1],[29,28],'k','LineWidth',1);
plot(2.1*[1.01,1],[29,28],'k','LineWidth',1);
ylabel('\sigma (deg)'); xlabel('L*');
set(gca,'Fontsize',18);
tx = text(0.13,57.3,'c'); set(tx,'FontSize',30);
set(gca,'XScale','log');
leg = legend('Exp 1','Exp 2','Exp 3','Exp 4','Exp 5',...
    'fontsize',16);
set(gca,'Position',[0.57 0.15 0.38 0.80]);
% set(leg,'Position',[0.81 0.70 0.08 0.21]);


%% Saving

if saveArg == 1
    fileName = 'Distinguish erupt'; 
    print(gcf,fileName,'-djpeg','-r300');
    saveas(gcf,fileName,'fig');
    saveas(gcf,fileName,'epsc');
end




